#include <iostream>
#include <stdio.h>

using namespace std;

//typedef float f32_t;
//
//class Complex {
//public:
//	f32_t real, imag;
//
//	friend std::ostream& operator<<(std::ostream &s, Complex z) {
//		return s << '(' << z.real << ',' << z.imag << ')';
//	}
//	Complex() {}
//	Complex(const f32_t &r, const f32_t &i) : real(r), imag(i) {}
//	friend Complex operator+(const Complex &l, const Complex &r) {
//		return Complex(l.real + r.real, l.imag + r.imag);
//	}
//	friend Complex operator-(const Complex &l, const Complex &r) {
//		return Complex(l.real - r.real, l.imag - r.imag);
//	}
//	friend Complex operator*(const Complex &l, const Complex &r) {
//		return Complex(l.real * r.real - l.imag * r.imag,
//		               l.real * r.imag + l.imag * r.real);
//	}
//	friend Complex fma(const Complex &x, const Complex &y, const Complex &d) {
//		f32_t r, i;
//		r = (x.real * y.real) + d.real;
//		i = (x.real * y.imag) + d.imag;
//		r = -(x.imag * y.imag) + r;
//		i = (x.imag * y.real) + i;
//		return Complex(r, i);
//	}
//};
//
//class SU3 {
//public:
//	Complex e00, e01, e02,
//	        e10, e11, e12,
//	        e20, e21, e22;
//	SU3() {}
//	SU3(Complex e00, Complex e01, Complex e02,
//	    Complex e10, Complex e11, Complex e12,
//	    Complex e20, Complex e21, Complex e22) {
//		this->e00 = e00;
//		this->e01 = e01;
//		this->e02 = e02;
//		this->e10 = e10;
//		this->e11 = e11;
//		this->e12 = e12;
//		this->e20 = e20;
//		this->e21 = e21;
//		this->e22 = e22;
//	}
//
//	friend SU3 operator+(const SU3 &x, const SU3 &y) {
//		return SU3(
//		           x.e00 + y.e00, x.e01 + y.e01, x.e02 + y.e02,
//		           x.e10 + y.e10, x.e11 + y.e11, x.e12 + y.e12,
//		           x.e20 + y.e20, x.e21 + y.e21, x.e22 + y.e22
//		       );
//	}
//
//	friend SU3 operator-(const SU3 &x, const SU3 &y) {
//		return SU3(
//		           x.e00 - y.e00, x.e01 - y.e01, x.e02 - y.e02,
//		           x.e10 - y.e10, x.e11 - y.e11, x.e12 - y.e12,
//		           x.e20 - y.e20, x.e21 - y.e21, x.e22 - y.e22
//		       );
//	}
//	
//	friend SU3 operator*(const SU3 &x, const SU3 &y) {
//		Complex t00, t01, t02,
//			t10, t11, t12,
//			t20, t21, t22;
//		t00 = fma(x.e00, y.e00, fma(x.e01, (y.e10), x.e02 * (y.e20)));
//		t01 = fma(x.e00, y.e01, fma(x.e01, (y.e11), x.e02 * (y.e21)));
//		t02 = fma(x.e00, y.e02, fma(x.e01, (y.e12), x.e02 * (y.e22)));
//		t10 = fma(x.e10, y.e00, fma(x.e11, (y.e10), x.e12 * (y.e20)));
//		t11 = fma(x.e10, y.e01, fma(x.e11, (y.e11), x.e12 * (y.e21)));
//		t12 = fma(x.e10, y.e02, fma(x.e11, (y.e12), x.e12 * (y.e22)));
//		t20 = fma(x.e20, y.e00, fma(x.e21, (y.e10), x.e22 * (y.e20)));
//		t21 = fma(x.e20, y.e01, fma(x.e21, (y.e11), x.e22 * (y.e21)));
//		t22 = fma(x.e20, y.e02, fma(x.e21, (y.e12), x.e22 * (y.e22)));
//		
//		return SU3(t00, t01, t02,
//			t10, t11, t12,
//			t20, t21, t22);
//	}
//};
//
//class ComplexMat3x3
//{
//public:
//	Complex _elems[3][3];
//	
//};
//
//int main()
//{
//	Complex a(1.0, 1.0), b(2.0, 2.0);
//	cout << a + b << endl;
//	cout << a - b << endl;
//	cout << a * b << endl;
//	cout << fma(a, b, b) << endl;
//	
//	return 0;
//}


template<class T, int N=3, int M=3>
struct mat {
	union {
		T a[N*M];
		T aa[N][M];
	}m;
	
	friend void mul(mat &result, const mat &a, const mat &b) {
		for(int i=0; i<N; i++) {
			for(int j=0; j<M; j++) {
				result.m.aa[i][j] = 0;
				for(int k=0; k<N; k++) {
					result.m.aa[i][j] += a.m.aa[i][k] * b.m.aa[k][j];
				}
			}
		}
	}
};

template<class T, int N, int M, int K>
void matmul(mat<T, N, M> &result, const mat<T, N, K> &a, const mat<T, K, M> &b) {
	for(int i=0; i<N; i++) {
		for(int j=0; j<M; j++) {
			result.m.aa[i][j] = 0;
			for(int k=0; k<K; k++) {
				result.m.aa[i][j] += a.m.aa[i][k] * b.m.aa[k][j];
			}
		}
	}
}


int main()
{	
	const int N = 16;
	const int M = 16;
	const int K = 16;
	
	mat<float, N, K> a;
	mat<float, K, M> b;
	mat<float, N, M> c;
	
//	for(int i=1; i<=N*K; i++) {
//		a.m.a[i-1] = i;
//		printf("%8d%s", i, (i%K==0 ? "\n" : "  "));
//	}
	for(int i=1; i<=N; i++) {
		for(int j=1; j<=K; j++) {
			a.m.aa[j-1][i-1] = (i-1)*N+j;
		}
	}
	for(int i=1; i<=N*K; i++) {
		printf("%.f%s", a.m.a[i-1], (i%K==0 ? "\n" : "  "));
	}
	cout << endl;
	for(int i=1; i<=K*M; i++) {
		b.m.a[i-1] = i;
		printf("%8d%s", i, (i%M==0 ? "\n" : "  "));
	}
	for(int i=1; i<=N*M; i++) {
		c.m.a[i-1] = 0;
	}
	cout << endl;
	matmul(c, a, b);
	
	for(int i=0; i<M; i++) {
		for(int j=0; j<N; j++) {
			printf("%.f  ", c.m.aa[i][j]);
		}
		cout << endl;
	}
	
	return 0;
}
